Geostatistics is a discipline concerned with spatially distributed random variables (also called “regionalized variables”), usually applied to problems in the earth sciences, such as estimation of mineral reserves and delineation of mineral deposits, hydrocarbon reservoirs, and groundwater aquifers. Typically, geostatistics makes use of two-point statistics summarized in a variogram. Multipoint (or multiple-point) geostatistics (MPS) differs from other variogram-based geostatistics primarily in that it characterizes spatial variability using patterns (sets of points and their configurations) that involve higher order (much greater than order 2) statistics.
One of the goals of multipoint geostatistics is simulation, namely the generation of numerical values along a line, on a surface, or in a volume, such that the set of values match certain given spatial correlation or structural properties (usually derived from a data set called an “analog” or “training image” or “template”) while optionally (in the case called “conditional simulation”) matching predetermined data. In practice, the “analog” may be, for example, a well-known rock volume that is statistically similar to a currently uncharacterized oil reservoir being delineated, and the predetermined data to be matched may be lithology observations at wells, or probabilities of lithologies derived from seismic data. MPS simulations are developed to match two types of constraints: absolute constraints are matched exactly; partial constraints are matched probabilistically, as closely as possible to the constraint value, but they are not guaranteed to match exactly. In the case of petroleum reservoir modeling, examples of absolute constraints are typically data acquired in wells or geological outcrops. Partial constraints are typically derived from seismic data; 1D, 2D, or 3D interpreted spatial trend information; geological category probability fields; or rotation and affinity (or scale) constraints on the size and direction of geological features. Such data are used in a stochastic modeling process to generate one-dimensional (ID), two-dimensional (2D), and/or three-dimensional (3D) spatial distribution of geological categories or rock properties. Since there is a random component involved in MPS simulations, individual stochastic realizations of property fields created by MPS algorithms differ, but the ensemble of realizations provide modelers with improved quantitative estimates of the spatial distribution and uncertainty of values in a modeled volume of interest.
Multipoint geostatistical methods often use a numerical training image to represent the spatial variability of geological information. The training image provides a conceptual quantitative description of the subsurface geological heterogeneity, containing possibly complex multipoint patterns of geological heterogeneity. Multipoint statistics conditional simulation anchors these patterns to well data (and/or outcrop data) and to the seismic-derived information (and/or probability field information or constraint grid(s)).
One multipoint geostatistics method, referred to as “Single Normal Equation Simulation” (SNESIM), is generally considered useful for practical applications such as modeling categorical or discrete data types, especially for categorical data in 3D property modeling. In the SNESIM method, the conditional probability density function of all categories at one point is computed using knowledge of the value at a number of nearby points and statistics provided by the training image. SNESIM works with discrete values only (i.e., a finite and usually small number of categories, such as for example five different rock types). Assume there are two categories to be simulated: A (“non-channel”) and B (“channel”). The training image contains a complete representation (i.e., an example) of a spatial distribution of A and B. Assume further that the category present at 4 points, which could represent wells, is known. Conceptually, the SNESIM method computes the probability of finding categories A or B at an unsampled point by scanning the training image for all occurrences of the “pattern” (that is, the pattern is defined by the relative spatial locations of all five points, and by the categories found at the four points where a value already has been assigned). If five such pattern occurrences (called replicates) are found, and 4 out of 5 replicates show category B at the relative position of the unknown point, then the method computes the probability of B at the unknown point, given the particular five-point pattern, to be ⅘ or 80% (while that of A is set to ⅕ or 20%). Furthermore, the method can assign a category to the unknown point by randomly drawing a value from a distribution with 80% probability of B and 20% probability of A if there are no replicates of the four-point pattern or its sub-patterns found in the training image.
In practice, the SNESIM uses pixel-based (or voxel-based in 3D) sequential simulation method. It starts with a volume to be modeled where all property values are unassigned at all grid cell locations, or one that contains only a few data points to be matched. These volumes are usually represented by a Cartesian grid, where each grid location is called a cell. First, SNESIM decides on a random path for visiting each unassigned cell once and only once. In the first cell, the method searches for nearby known points within a search volume, usually an ellipsoid or rectangular volume around the known point. If one or more known (or already assigned) cells are found, it proceeds as to the next point on the random path as described above to find the probability of finding categories A or B at the unknown point. Armed with the probabilities for each category at the unknown point, the method randomly draws a value (category A or B) from the inferred probability distribution and assigns it to the unknown point. The process repeats at the next cell in the initially assigned random path until all cells are assigned.
The SNESIM method also provides an advantage over other MPS algorithms in that it reduces computations by storing the information for a training image in a tree-like structure to help search for a pattern in the training image. The pattern is only known on a subset of points within a neighborhood of the points to be simulated, which means that the set of locations in the training image that match the pattern generally cannot generally be found with one single search in the tree.
One alternative approach to SNESIM, referred to as a Markov mesh method, is to use a structured path instead of a random path for visiting unassigned cells. Similar to SNESIM, Markov mesh methods emulate one dimensional stationary Markov chains in that, when simulating a point, they have a ‘history’ of points which have been simulated already, and a ‘future’, which are the points that have not yet been simulated. The search neighborhood of a point to be simulated selects a subset of the points from the history. The major differences with SNESIM are due to the fact that the set of points in the search neighborhood are the same for each point to be simulated, which is important in two ways. Firstly, this makes the random function stationary, which helps with the reproduction of the training image patterns in that smaller search neighborhoods can often be used, thereby making the simulation faster. Secondly, the search neighborhood is full, which makes the tree search faster than the case of a random path, where the search needs to account for the missing cells in the search neighborhood. A unilateral path method ensures that only one traversal of the tree is needed to find the examples of the pattern within the training image as opposed to the multiple traversals needed with the random path method. The most significant problem for unilateral methods is that it is far more difficult to condition to well data, as well data points can occur in the future of points that are to simulated. While this may not be a significant problem for typical applications of unilateral methods in image analysis that do not have to account for data, it is a significant issue for applications to geological modeling where conditioning to data, in the form of wells or observations, is important.
Furthermore, the SNESIM method visits each unassigned cell in a random path, which in practice results in search neighborhoods, or a “search mask” with many unassigned cells. When an unassigned cell lies within a search mask, the calculation of the conditional probability becomes more computationally intensive due to traversal of more paths within the tree.
Moreover, the SNESIM method shows a computational drawback when simulating cells where the search neighborhood extends past the simulation grid, as in practice the method assumes all cells outside of the grid are equivalent to unassigned, not-yet-simulated, cells. Again the introduction of unassigned cells into the search neighborhood dampers the performance of the algorithm.
The SNESIM method is also typically limited to single-threaded execution. While computer technology has advanced in the area of parallelization, such that many computer workloads can be broken down into component pieces and executed in parallel using multiple hardware threads executing in multiple processing cores, processors, processing nodes and/or networked computers, the SNESIM method is not readily adapted to parallelization due to the interdependence of cell simulations. Given that the assignment of a cell is dependent upon the values assigned to neighboring cells, substantial synchronization of cell simulations would be required to ensure that individual cell simulations were based on accurate simulations of neighboring cells, effectively leading to substantial delays in threads while waiting on the results from other threads, and negating any potential gains from devoting multiple hardware threads to the simulation.
Therefore, a substantial need continues to exist in the art for an improved manner of performing geostatistics models that is less computationally intensive than conventional approaches.